Signature of Nearly Icosahedral Structures in Liquid and Supercooled Liquid Copper 
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A growing body of experiments display indirect evidence of icosahedral structures in supercooled 
liquid metals. Computer simulations provide more direct evidence but generally rely on approximate 
interatomic potentials of unproven accuracy. We use first-principles molecular dynamics simulations 
to generate realistic atomic configurations, providing structural detail not directly available from 
experiment, based on interatomic forces that are more reliable than conventional simulations. We 
analyze liquid copper, for which recent experimental results are available for comparison, to quantify 
the degree of local icosahedral and polytetrahedral order. 

PACS numbers: 61.43.Dq,61.20. Ja,61.25.Mv 
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I. INTRODUCTION 

TurnbuU |1{ established that metallic liquids can 
be supercooled if heterogeneous nuclcation can be re- 
duced or avoided. Later, Frank hypothesized that the 
supercooling of liquid metals might be due to frustrated 
packing of icosahedral clusters. Icosahedral clustering 
of 12 atoms about a sphere is energetically preferred to 
crystalline (e.g. FCC, HCP or BCC) packings for the 
Lennard- Jones (L-J) pair potentials. The icosahedron 
is favorable because it is made up entirely of four-atom 
tetrahedra, the densest-packed cluster possible. Local 
icosahedral order cannot be propagated throughout space 
without introducing defects. 

Remarkably, the frustration of packing icosahedra is re- 
lieved in a curved space, where aperfect 12-coordinated 
icosahedral packing exists |1, H, Q . Disclination line de- 
fects must be introduced into this icosahedral crystal in 
order to "flatten" the structure and embed it in ordinary 
three dimensional space. Owing to the 5-fold rotational 
symmetry of an icosahedron, the disclination lines are of 
type 72°. The negative disclination line defects that are 
needed to flatten the structure cause increased coordina- 
tion numbers of 14, 15 or 16. Large atoms, if present, 
would naturally assume high coordination number and 
aid in the formation of a disclination line network. 

Many studies of Lennard- Jones systems have tested 
Frank's hypothesis. Hoare found that for clusters 
ranging between 2 to 64 atoms at least three types of 
"polytetrahedral" noncrystalline structures exist, with a 
higher binding energy than HCP or FCC structures with 
the same number of atoms. Honeycutt and Andersen gj 
found the crossover cluster size between icosahedral and 
crystallographic ordering around a cluster size of 5000 
atoms. They also introduced a method to count the num- 
ber of tetrahedra surrounding an interatomic bond. This 
number is 5 for local icosahedral order. Steinhardt, Nel- 
son and Ronchetti introduced the orientational order 
parameter VPg to demonstrate short range icosahedral or- 
der. 

Many other simulations have been performed on pure 



elemental metals and metal alloys, using a modified John- 
son potential [^, embedded atom potenti als ll, .12j. the 
Sutton-Chen (SC) many body potential [13, to name 
a few. These potentials model the interatomic interac- 
tions with varying, and generally uncontrolled, degrees 
of accuracy. Ab-initio studies on liquid Copper [l4.[l5|. 
Aluminum and Iron ^3 achieve a high degree of re- 
alism and accuracy, but have not been analyzed from the 
perspective of icosahedral ordering. Nevertheless, recent 
Ab-initio studies on Ni and Zr |0, ^ have been done 
with this perspective and find that with supercooling the 
degree of icosahedral ordering increases in Ni while in Zr 
BCC is more favored. 

X-ray diffraction measurements of electrostatically lev- 
itated droplets of Ni "2^ found evidence of distorted 
icosahedral short ranged order. Neutron scattering stud- 
ies of deeply undercooled metallic melts |21| observed 
the characteristic shoulder on the second peak of the 
structure factor, which has been identified as a signature 
of icosahedral short range order |2^ |2^ . The shoulder 
height increases with decrease in temperature. 

A recent XAS (X-ray Absorption Spectroscopy) ex- 
periment on liquid and undercooled liquid Cu by Di Ci- 
cco et al. |2^ isolated the higher order correlation func- 
tions. They a ppl ied Reverse Monte-Carlo (RMC) refine- 
ment |2^ |2y, [23 simultaneously to diffraction and XAS 
data to construct a model of the disordered system com- 
patible with their experimental data. They analyze the 
three body angular distribution function N(0) and also 
the orientational order parameter W^g- Their conclusion 
was that weak local icosahedral order could be observed 
in their sample. This experiment provided the most di- 
rect experimental evidence to-date of the existence of 
icosahedra in a liquid metal. 

Motivated by these results, we explore the structures 
of liquid and undercooled liquid metals using first prin- 
ciples simulations. First principles calculations achieve 
the most realistic possible structures, unhindered by the 
intrinsic inaccuracy of phenomenological potentials, and 
with the ability to accurately capture the chemical na- 
ture and distinctions between different elements and al- 
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loys. We use the VASP (Vienna Ab-initio Simulation 
Package code which solves the quantum me- 

chanical interacting many-body problem using electronic 
density functional theory. These forces are incorporated 
into a molecular dynamics simulation. The trade-off for 
increased accuracy is a decrease in the system sizes we 
can study, so we can only observe local order, not long 
range. Also we are limited to short time scales. 

Our analysis covers methods that have previously been 
fruitful. We look at the radial distribution function, 
the structure factor, the three body angular distribution 
function, which is simply related to the three body cor- 
relation function, the Wq parameter as discussed above, 
and the Honeycutt and Andersen bond statistics method 
i 

The extent of the icosahedral order that we observe 
in simulation is qualitatively in agreement with recent 
experiments (2^. At high temperatures we found that 
structural properties of liquid Cu strongly resembled a 
maximally random jammed |3C| hard sphere configura- 
tion. From this we conclude that a nearly universal struc- 
ture exists for single component systems whose energet- 
ics are dominated by repulsive central forces. The degree 
of icosahedral order is not great, presumably due to the 
frustration of icosahedra, but it does show a tendency to 
increase as temperature drops. 

Section ^] describes our first principles molecular dy- 
namics method in greater detail. The next section, sec- 
tion lllll discusses our study on copper. Here we introduce 
the radial distribution function g(r), the liquid structure 
factor S{q), the We bond orientational order parameter, 
the three body angular distribution function N{0), and 
the Honeycutt and Andersen analysis method. We con- 
clude, in section llVl with some thoughts about enhancing 
icosahedral order by alloying with a fraction of smaller 
and larger atoms. 



II. FIRST PRINCIPLES METHOD 

First principles simulation is an incisive, powerful 
and well-developed tool based on a quantum mechani- 
cal treatment of the electrons responsible for interatomic 
bonding. Since the method is based on fundamental 
physical laws and properties of atoms, it can be applied 
to a wide variety of metals, including alloys, and yields 
the energy and forces computationally without any ad- 
justable free parameters. 

Our ab — initio molecular dynamics simulation pro- 
gram, VASP 28^ solves the N-body quantum me- 
chanical interacting electron problem using electronic 
density functional theory, under the Generalized Gra- 
dient Approximation (GG A). We use the projector- 
augmented wave m (PAW) potentials as provided 
with VASP. Calculation times grow nearly as the third 
power of the number of atoms, limiting our studies to 
sample sizes of around a hundred atoms. 

In first-principles molecular dynamics, although inter- 



atomic forces and energies are calculated quantum me- 
chanically, we still treat the atomic motions classically, 
using the Bor n-Op penheimer approximation. We use 
Nose dynamics '3s|| to simulate in the canonical ensemble 
at fixed mean temperature. The system was well equili- 
brated before data was considered for analysis. The sim- 
ulation started with a random configuration, at a tem- 
perature high enough to ensure a liquid state, and was 
allowed to equilibrate at this high temperature. Subse- 
quently, lower temperatures were simulated starting from 
previous configurations. All calculations were F point 
calculations (a single 'k' point). 



We took N=100 Cu atoms and applied periodic bound- 
ary conditions in an orthorhombic cell. Our unequal lat- 
tice parameters avoid imposing a characteristic length on 
the system. The simulations were done at three different 
temperatures, T=1623K, 1398K and 1313K in order to 
compare with Di Cicco's experiments (24] . The melting 
point of copper is T=1356K, so samples at 1623K and 
1398K are in the liquid regime, while the one at 1313K is 
undercooled. We used number densities of 0.0740 A^'^, 
0.0758 A-3 and 0.0764 A'^ respectively at T=1623K, 
1398K and 1313K. These were obtained from a fit of the 
XRD experimental volume per particle f33| to a straight 
line versus temperature. Starting from configurations 
that had been previously equilibrated at slightly differ- 
ent densities, transients of about 250 steps (Ifs per step) 
passed prior to the onset of equilibrium fluctuations of 
the energy. After the transient, a total of 3000 MD steps 
were taken at each temperature, for a total simulation 
time of 3ps. The run time was around 480hrs on a 2.8 
GHz Intel Xeon processor for each temperature. Subse- 
quently, two configurations from the highest temperature 
run, widely separated in time, were selected and used as 
the starting configuration for two runs at T=1398K, with 
proper scaling of densities. After an initial transient of 
250 steps, these runs were further continued at T=1313K 
as weU as at T=1398K. The four runs, two at T=1398K 
and two at T=1313K were carried on for Ips, during 
which time the energy showed equilibrium fluctuations. 
All of these runs have been used for analyzing the local 
order in liquid copper. 



For a few selected configurations, a conjugate-gradient 
algorithm was used to relax the ions to their instan- 
taneo us g round state, to explore their inherent struc- 
tures [35|- Surprisingly all of our samples partially crys- 
tallized, based on visual observation. Further efforts 
to obtain quenched amorphous structures used steepest- 
descent minimization and molecular-dynamics with a lin- 
ear temperature ramp, followed by steepest-descent min- 
imization. Again the structures partially crystallized. 
There appear to be no fully amorphous relaxed struc- 
tures accessible for Cu using these methods. 
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III. RESULTS 

A. Radial Distribution Function g{r) 

The radial distribution function, g{r), is proportional 
to the density of atoms at a distance r from another 
atom. We calculate g{r) by forming a histogram of bond 
lengths. We use the repeated image method to obtain 
bond lengths greater than half the box size, and antic- 
ipate g(r) in this range may be strongly influenced by 
finite size effects. We then smooth out with a gaussian 
of standard deviation 0.05 A. Fig. ^ shows the g(r) we 
obtained at the three different temperatures. Our g(r) 
at T=1623K, compares well with g{r) interpolated from 
XRD experiments 0|, with the two curves overlapping 
almost everywhere except for a small disagreement in the 
position of the first peak. Results from neutron diffrac- 
tion experiment [sBl at T=f393K, compare well with our 
g{r) at T=f398K. Comparisons with the g{r) for Cu at 
fSOOK from the ab — initio MD studies by Hafncr et. 
al. 01 and Vanderbilt et. al. 0| finds that the heights of 
their first peak match well with our g{r) (interpolated to 
T=1500K). But their peak positions are shifted slightly 
to the left of ours (ours is at 2.50 A). The g(r) from an 
embedded-atom method (EAM) model for Cu ,1^] which 
matches almost exactly with the XRD data at T=1773K, 
is also consistent with our extrapolated g{r) at this tem- 
perature. 

The growth in height of the peaks in the supercooled 
system at T=1313K suggests an increase of some type of 
order. However this order is not related to the crystalline 
FCC equilibrium phase, as we show in the following sub- 
sections. 

To test for finite size effects in our N=100 atom sys- 
tem, we ran a separate simulation for N=200 atoms at 
the intermediate temperature T=1398K (Fig. 13). The 
first and second peaks of g{r) for both the system sizes 
compare very well. There is a small but significant differ- 
ence in the depth of the first minimum, then systematic 
differences between the curves beyond 5A. From this we 
conclude that the finite size effect is not important at 
small r values, but for larger values of r (beyond 5 A) 
there is a weak finite size effect. The three body angular 
distributions, and the Wq histograms of the N=200 and 
the N=fOO runs, are also comparable, suggesting that 
N=100 is sufficient for studies of local order of the types 
we consider here. 

We calculate the coordination number from the radial 
distribution function g{r). We choose a cutoff distance 
near the first minimum of g(r), at i?cMt=3.4 A. The pre- 
cise location of the minimum is difficult, and its variation 
with temperature is smaller than the error in locating its 
position (Fig. so that we don't change the value of 
Rcut with temperature. With this value of Rcut we find 
an average coordination number (iVc) of f2.3 which is 
nearly independent of temperature {Nc changes from f 2.f 
at high temperature to 12.5 with supercooling). However, 
our range of evolution of Nc is small compared to the case 
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FIG. 1: (Color online) Simulated liquid Cu radial distribu- 
tion function, g{r), at three different temperatures. The 
simulated curve at T=I623K matches well with the exper- 
imental XRD (X-Ray Diffraction) result [s^ interpolated to 
T=1623K. (The simulated g{r) at T=f3f3K and T=I398K 
have been shifted up for visual clarity) 

of Ni |l8l | , since our degree of supercooling is much lower 
(3%) than theirs (17%). We are not able to achieve a 
higher degree of supercooling of Cu as mentioned earlier 
in the paper. 

B. Liquid Structure Factor S(q) 

The liquid structure factor S{q) is related to the radial 
distribution function g{r) of a liquid with density p by, 

oo 

S{q) = l + 4Trp f [g{r) - Ij^^^^r^dr. (1) 

J ' qr 



Evidently, one needs a knowledge of the radial distribu- 
tion function up to large values of r to get a good S{q). In 
our first principles simulation, we are restricted to small 
values of r, due to our small system sizes, so we need a 
method to get S{q) from our hmited range g{r) function. 

Baxter developed a method 'st', 'ss'l to extend g{r) be- 
yond the size of the simulation cell. The method exploits 
the short ranged nature of the direct correlation func- 
tion c(r), which has a range similar to the interatomic 
interactions 39], as opposed to the g{r) which is much 
long ranged. The exact relation that connects these two 
fmictions is the Ornstein-Zernike relation, 

h{r) = c(r) +p j h{\v- r'|)c(|r'|)dr' (2) 

where h{r) — g{r) — 1. 
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FIG. 2: (Color online) Simulated liquid Cu radial distribution 
function, g{r), at T=1398K for N=100 and N=200 atoms. A 
weak finite size effect is observed after about r—sA. 



Assuming that c(r) vanishes beyond a certain cutoff 
distance rc, Baxter obtained a pair of equations, vahd 
for r < re- The remarkable property of this method is 
that if we know h{r) over a range < r < Vc, then we 
can obtain c(r) over its entire range (from to Tc), which 
imphcitly determines h{r) over its entire range (from 
to oo) through Eq. 121. 

We solve the Baxter's equations iteratively to obtain 
the full direct correlation function. A complete knowl- 
edge of the direct correlation function gives us the struc- 
ture factor S{q) in terms of its fourier transform c{q), 



S{q) 



1 - pc{q) 



where, 



qr 



(3) 



(4) 



The S{q) showed good convergence with different choices 
of Tc, and a choice of rc=5A seemed appropriate because 
it was one half of our smallest simulation cell edge length. 
Even though in metals there are long range oscillatory 
Friedel oscillations, our ability to truncate c(r) at rc—5A, 
shows that these are weak compared with short range 
interactions. 

Fig. Q compares the calculated S{q) at our three dif- 
ferent temperatures, and the experimental neutron S{q) 
at T=1393K The calculated S{q) at T=1398K com- 
pares well with the experiment at all values of q. The 
S{q) from the larger system is in better agreement with 
the experiment. Comparison between the two system 



FIG. 3: (Color online) Liquid structure factor S{q) as ob- 
tained from the simulated radial distribution function g{r) at 
T=1398K compared with the S{q) from neutron diffraction at 
T=1393K (3^ . The calculated S{q) at the other two temper- 
atures are also plotted, and show the expected temperature 
behavior. (The S{q) at T=1313K and T=1398K have been 
shifted up for visual clarity) 



sizes suggest again that the finite size effects are signif- 
icant but not important, and N=100 is good enough to 
get a representative liquid structure. No resolution cor- 
rection was applied to the experimental data, and more- 
over it was smoothed. Both of these cause a decrease in 
the height of the actual S{q), which becomes quite ap- 
preciable at the first peak. As a test, we also applied a 
resolution correction to our simulated S{q) (not shown), 
which reduced the height of the first peak bringing it 
in closer agreement with the experimental value. Nev- 
ertheless the overall excellent agreement shows that the 
first principles simulation with only N=100 atoms is able 
to produce representative structures at T=1398K. This 
enables us to make further studies of the local icosahe- 
dral and polytetrahedral order in liquid and supercooled 
liquid copper. 

As mentioned earlier in the introduction, one signature 
of icosahedral short range order is the splitting of the sec- 
ond peak of S{q) 22, 23|. Even though we observe weak 
icosahedral order in Cu (discussed later in this paper), 
we do not observe a clear splitting of the second peak in 
S{q) (Fig. El, but we do observe a broadening as we lower 
the temperature. We think that the absence of splitting 
could be because of our low degree of supercooling. 



C. Bond Orientation Order Pcirameters 

As introduced by Steinhardt, et al. 0, the Wi param- 
eters are a measure of the local orientational order in 
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TABLE I: We values for a few clusters 
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FIG. 4: (Color online) Simulated Wa distributions for liquid 
Cu. Ideal icosahedron and FCC values are indicated. 



liquids and undercooled liquids. To calculate Wi, the 
orientations of bonds from an atom to its neighboring 
atoms are projected onto a basis of spherical harmon- 
ics. Rotationally invariant combinations of coefficients in 
the spherical harmonics expansion are then averaged over 
many atoms in an ensemble of configurations. The re- 
sulting measures of local orientational order can be used 
as order parameters to characterize the liquid structures. 
For an ideal icosahedral cluster, ^ = 6 is the minimum 
value of I for which We ^ 0. Table enumerates Wq 
values for different ideal clusters. We see that the ideal 
icosahedral value of We is far from other clusters, making 
it a good icosahedral order indicator. 

We choose the cutoff distance to specify near neigh- 
bors at Rcut—SA A as before. Our value of Rcut is 
significantly greater than that of Di Cicco, et ai. Our 
Wq distributions (Fig. Q show strong asymmetry favor- 
ing negative values with tails extending towards the ideal 
icosahedron value. Because the histogram vanishes as Wq 
approaches its limiting negative value we see that there 
are essentially no perfectly symmetric undistorted icosa- 
hedra present in our simulation. However, a significant 
fraction do have Wq values close to the icosahedral value. 

A We analysis was performed for a 10^ atom maximally 
random jammed hard sphere configuration [s^- The di- 
ameter of the hard spheres was rescaled so that the posi- 
tion of the main peak of the resulting g{r) matched the 
value r = 2.5A found for Cu at T=1623K. The Rent 
value for the MRJ configuration was taken near the first 
minimum of the g{r) at r = 3.3A. Remarkably, the Wq 
distribution of the MRJ configuration (Fig. 0)) is simi- 
lar to the distribution for liquid Cu at high temperature, 
suggesting that the structure of Cu under this condition 
is dominated by strongly repulsive short-range central 
forces. 
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BCC 
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12 


12 


12 
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We 


-0.012 


-0.013 


-0.169 


+0.013 



As we lower the temperature of liquid Cu, the mean 
value of We drops and the peak of the We distribution 
shifts to the left. However, the peak never moves be- 
low We = -0.05, and the tail of the distribution at neg- 
ative We shows no strong temperature dependence. It 
seems that there is no change in the number of nearly 
icosahedral clusters as the temperature drops into the 
supercooled regime, possibly a result of the frustration 
of icosahedral packing. Our liquid has a single compo- 
nent, so there is no natural way to introduce disclina- 
tions. This inhibits the growth of a population of atoms 
with We close to its ideal icosahedral value. 

Comparing our result with that of Di Cicco et al. at 
T=1313K, we see that our curve is more asymmetric to- 
wards negative values than Di Cicco's, so that we see a 
greater fraction of atoms near the ideal icosahedral value 
of We- The discrepancy probably lies in the difference 
between the two methods used to generate the positional 
configurations (the difference is even greater if we use 
Di Cicco's value of Rcut)- Their configurations were ob- 
tained using Reverse Monte Carlo (RMC), which does 
not guarantee accurate configurations. Our first prin- 
ciples method should be more accurate in determining 
these configurations. Of course, Di Cicco's configurations 
are consistent with experimentally measured three-body 
correlations. It would be of great interest to see if our 
configurations are also consistent with the raw experi- 
mental data. The differences in We distributions should 
not be overstated - the experiment and our simulations 
both show that liquid and supercooled liquid copper has 
weak but non-negligible icosahedral order. 



D. Bond Angle Distribution N(6): 

The bond angle distribution N{6) is a simple type of 
three-body correlation function. Let be the angle be- 
tween bonds from a single atom to two neighbors, and 
define N(6') as the probability density for angle 9, normal- 
ized such that the total probability, J N{9)dd—1. The 
distribution for the central atom of an ideal 13-atom 
icosahedral cluster, shows peaks at 63.4°, 116.4° and 
180.0°. For other crystallographic clusters, like HCP, 
FCC, and BCC, we see peaks at 60°, 90° and 120°. An- 
gles around 60° degrees indicate nearly equilateral trian- 
gles that may well belong to tetrahedra. 

Fig.|31shows the distributions for copper at three tem- 
peratures. We have chosen the same value of Rcut that 
was used to obtain We in section lTlI CI . The distribution 
function shows maxima at 56° and 110° with a minimum 
around 80° . Our result is similar to that of Di Cicco (they 
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FIG. 5: (Color online) Distribution of N(6l) for liquid Cu. 
Ideal icosahedron values are indicated. 



show only T=1313K), but with more pronounced mini- 
mum and second maximum. The peak around 60° shows 
an abundance of nearly equilateral triangles, indicating 
the presence of tetrahedrons, which can pack to form 
icosahedra. The minimum close to 90° shows that there 
aren't many cubic clusters. We also see that the high- 
angle tail at high temperature turns into a broad max- 
imum at low temperature centered around 165°. This 
may represent a shifting of the ideal 180° peak caused by 
cluster distortion. The ordering increases as temperature 
decreases, indicating that the number of nearly equilat- 
eral triangles increases when the liquid is undercooled, 
probably caused by an increase in polytetrahedral order 
with undercooling. 

The distribution of the MRJ configurations (same Rcut 
as defined in section IIII ("P shows a sharp peak at ex- 
actly 60°, a broad peak at 110° and a minimum around 
90°. The peak at 60° shows an overwhelming presence 
of perfectly equilateral triangular faces, which are easily 
formed when 3 hard spheres come in contact with each 
other. But the minimum around 90° and a second max- 
imum nearer to 110° as opposed to 120°, suggests that 
the local order is not FCC or HCP. This feature of the 
MRJ configuration agrees qualitatively with the angular 
distribution of liquids, implying an underlying universal 
structure for systems whose energetics are dominated by 
repulsive central forces. But the quantitative differences 
also emphasize the necessity to exactly model an atomic 
liquid to study its local environments, and quantify poly- 
tetrahedral order. 



FIG. 6: (Color online) Honeycutt- Andersen pair fractions for 
142's (FCC and HCP forming +72° disclination) , 15's (icosa- 
hedron) and 16's (—72° disclination) at different tempera- 
tures. The increase in 15's as temperature is reduced show 
increased icosahedral ordering with supercooling. The corre- 
sponding HA values for the MRJ configurations are indicated 
on the right side of the plot. 



E. Honeycutt and Andersen analysis 

Honeycutt and Andersen [|| introduced a useful assess- 
ment of local structure surrounding interatomic bonds. 
We employ a simplified form of their analysis, count- 
ing the number of common neighbors shared by a pair 
of near-neighbor atoms. This identifies the number of 
atoms surrounding the near-neighbor bond and usually 
equals the number of edge-sharing tetrahedra whose com- 
mon edge is the near-neighbor bond. We assign a set of 
three indices to each bond. The first index is 1 if the 
root pairs are bonded (separation less than or equal to 
Rcut)- The second index is the number of near-neighbor 
atoms common to the root pairs, and the third index 
gives the number of near-neighbor bonds between these 
common neighbors. We take i?cMt=3.4A as before. Note 
that the Honeycutt and Andersen fractions depend sensi- 
tively on Rcut, making precise quantitative comparisons 
with other prior studies difficult. 

In general, 142's are characteristic of close packed 
structures (FCC and HCP) and 143's are characteristic 
of distorted icosahedras ;4QJ. They can also be consid- 
ered as +72° disclinations 0, 0, |g. Likewise, 15's are 
characteristic of icosahedra, and 16's indicate -72° discli- 
nations. Fig.|H|shows the 14's, 15's and the 16's for liquid 
Cu at the three temperatures. The error bars shown were 
calculated by breaking the data into three subsets. The 
14's have been separated into 142's and 143's. The re- 
maining 14's are mostly 144's with fraction around 0.04. 
The fraction of 142's and 143's holds steady with tem- 
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perature, with the icosahedral fraction always exceeding 
the close packed fraction. As the temperature drops, the 
fraction of 15's grows. At each of the three temperatures, 
the 15's are mainly comprised of 154's ( characteristic of 
distorted icosahedra) and 155's (characteristic of perfect 
icosahedra), with the 154's slightly higher than the 155's. 
Of all the 16's, the 166's are the highest and steadily in- 
crease with lowering of temperature. The 166's indicate 
the -72° disclination lines, which relieve the frustration 
of icosahedral order. 

These trends indicate a weak increase in polytetrahe- 
dral ordering with supercooling. The same trend was 
observed in simulations based on Sutton-Chen poten- 
tials except for the fact that our 142's are slightly 
higher compared to their 142's. 

For comparison, the values for a maximally random 
jammed packing jSfl (i?cMt=3.3A) are shown in Fig. 
These values are fairly close to liquid Cu at high temper- 
ature, and also to a similar common neighbor analysis 
of dense random-packing of hard spheres [42| . These re- 
sults are consistent with our previous observation for the 
We distribution and N(0), that a nearly universal struc- 
ture arises at high temperature, dominated by repulsive 
central forces. 



IV. CONCLUSION 

This study quantifies icosahedral and polytetrahedral 
order in supercooled liquid copper. While the structural 



properties of high temperature liquid Cu are close to a 
maximally random jammed structure (3^], proper mod- 
eling of atomic interactions is essential to capture the 
behavior of an element at liquid and supercooled temper- 
atures. A first-principles simulation is the most reliable 
means of achieving this. We find small but significant dis- 
agreement with analysis based on Reverse Monte-Carlo 
simulation. 

Supercooled liquid copper shows a slight increase in 
icosahedral and polytetrahedral order as temperature 
drops, which is consistent with recent experiments poL 
I21L |24 |. The frustration of icosahedrons in the one 
component liquid inhibits formation of perfect icosahe- 
dra, giving rise to defective icosahedrons. Alloying with 
larger atoms might relieve the frustration of packing 
icosahedrons by encouraging the formation of -72° discli- 
nations Alloying with smaller atoms can re- 

lieve frustration of individual icosahedrons by placing the 
smaller atom at the center 0. Alloying with larger and 
smaller atoms simultaneously thus offers the chance to 
optimize icosahedral order. Work is in progress in achiev- 
ing the same l45l| . 
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